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The ensemble of now more than 250 discovered planetary systems displays a 

a . 

wide range of masses, orbits and in multiple systems, dynamical interactions. 
These represent the endpoint of a complex sequence of events, wherein an en- 
^ | tire protostellar disk converts itself into a small number of planetary bodies. 

Here we present self-consistent numerical simulations of this process, which 

O 1 

produce results in agreement with some of the key trends observed in the prop- 
CD ■ 

j>! | erties of the exoplanets. Analogues to our own Solar System do not appear to 

be common, originating from disks near the boundary between barren and 



(giant) planet-forming. 

Gas giant planets form within 1-10 million years (Myrs), during the time that their par- 
ent star possesses a gas disk (1). Two- and three-dimensional hydrodynamic simulations of 
planets embedded in protostellar disks are, due to their computational cost, limited to follow- 
ing a maximum of 10 3 — 10 4 orbits; even for a planet as far out as Jupiter (orbital period 11 
years), this amounts to at most a tenth of a protostellar disk's total lifetime. For this reason, 
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longer-timescale simulations usually model the time after the protostellar gas disk has dissi- 
pated, allowing the problem to be tackled with pure gravitational iV-body methods. Recent 
work in this area has been successful at reproducing the observed exoplanet eccentricity distri- 
bution, the key requirement being simply that planets must begin close enough to each other 
to render systems dynamically active from the outset (|2][£|). Although these results provide 
support for planet-planet interaction as the main agent behind the eccentricities, they do not 
address the question of how, or even if, planet formation produces the requisite crowded sys- 
tems. Furthermore, many of the observed cases appear to have kept a memory of their early 
dynamical evolution: About a quarter of discovered systems contain planet pairs locked into 
dynamical mean-motion resonances (@), likely the product of early differential migration © 
via planet-disk interaction (00©. Thus, to achieve a more complete understanding of the 
planet formation process, it is essential to bridge the disk and post-disk era. Most models of this 
regime have been either semi-analytic treatments, which do not consider planet-planet interac- 
tions dPlESIZZlIZ!]), or iV-body simulations with parameterizations of disk effects ((721[7?|). Here, 
we perform simulations using a hybrid numerical scheme (175)) . that combines the N-body code 
SyMBA (Symplectic Massive Body Algorithm) ({ZD with a one-dimensional disk model that 
self-consistently interacts with, and accretes onto, the embedded planets. We began our sim- 
ulations with the appearance of the first small protoplanets (1CT 3 Me where Me is the Earth's 
mass) in a gas disk and ran them for up to half a billion years (Gyrs). Gas giants are assumed to 
form via core accretion (see Supporting Material for a more detailed description of the code). 

The fundamental question we sought to address is how the properties of a mature planetary 
system map to those of its birth disk. To this end, we performed simulations covering a range 
of disk parameters; Fig. 1 depicts one representative example, which illustrates some of the key 
effects at work in the planet formation process (see Supporting Material for an animated version, 
as well as additional examples). In this simulation, protoplanet growth initially proceeded in 
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an orderly manner; the dissipational effect of the gas disk dominated over the bodies' mutual 
perturbations, so that all orbits were kept nearly circular. Simultaneously, interaction with the 
disk also caused protoplanets to undergo orbital decay, referred to as "type I migration" ©. 
This effect will in general deposit some protoplanets at the inner disk edge (or onto the central 
star) before they have a chance to become gas giants (|70l) . but since migration speed scales 
with migrator mass and gas disk density, those protoplanets that grow slowly enough relative 
to the disk's depletion time avoid this fate (1771172)) . All protoplanets in the simulation (Fig. 1) 
fell into the latter category (see Supporting Material for counterexamples), although the orbit 
of the innermost one, originating at 2 AU, shrank to 1/5 its original radius by 1 Myr. Around 
this time, the first gas giant formed, scattering many of the neighboring protoplanets to high 
eccentricities as its mass rapidly grew. Its perturbation on the surrounding gas also became 
significant, and a deep annular gap opened around its orbit. With gas unable to readily cross 
its orbit, the planet was locked into the disk and carried along by the latter' s accretion flow, 
in what is termed "type II migration" ©. At 1.2 Myrs, the gap become a hole as the inner 
disk material drained away faster than the planet migrated; however the hole rapidly closed in 
again, depositing the planet at the inner disk edge at 1.24 Myrs. By this time two more gas 
giants began to form. Each in turn opened a hole in the disk, in the process cutting off its inner 
neighbor from the gas, so that ultimately only the outermost planet was directly interacting 
with — and accreting from — the disk. The outer planet migrated toward the middle one, and 
at 1.66 Myrs, the two become locked into a 3:1 mean-motion resonance (18) with each other, 
mutually increasing their eccentricities as they moved inward together. Migration and further 
planet accretion stalled as the disk density became low; the disk disappeared altogether a little 
after 4 Myrs, leaving behind a system consisting of three gas giants — the outer two still in 
resonance — plus an outer Neptune-class planet. This system remained stable for the rest of the 
simulation, which ended after 0.5 Gyr. 



3 



We performed a set of 100 simulations over a range of parameters. The properties of pro- 
toplanetary disks are only weakly constrained; two key parameters are initial gas disk mass 
M disk , and disk viscosity u, which determines how rapidly accretion onto the central star re- 
moves the bulk of the gas [photoevaporation may remove the last of the gas (1791) 1. Observations 
combined with modelling (l20l) suggest that, roughly, 10~ 2 Ms un < Mdi S k < 10 _1 Msun and 
1CT 3 < a < 1CT 2 , using the common parameterization ([271) v = ac s H, with c s the gas sound 
speed and H the disk scale height. We used this range of parameters for our simulations; the 
resulting planetary systems (Fig. 2) were each the product of a complex interplay between 
planet-disk and planet-planet effects, thus individual outcomes were highly stochastic. Never- 
theless, clear trends with Mdi S k and a are visible: In one extreme, low disk mass combined with 
high viscosity resulted in systems that produced no gas giants at all. In the other extreme, high 
disk mass combined with low viscosity results in the production of numerous gas giants; most 
underwent significant inward migration, and many acquired large eccentricities. 

We can understand these results in terms of two fundamental timescales in a planet-forming 
disk: One is the gas disk depletion time Tdisk, the other is the time to form the first gas giant, 
Tgiant- As shown in Fig. 3, we expect an initial burst of planet formation that spreads out from 
a particular radius, typically comparable to the Jupiter-Saturn region of our own Solar System 
(5-10 AU). As time passes, the delay between the birth of successive giants becomes longer and 
longer. Thus, the crucial factor determining how a given system's formation history will play 
out is the time during the gas disk's lifetime that this burst occurs. In cases with r giant > r disk 
(Fig. 3, lower right), the gas is removed before any gas giant has a chance to form, leaving 
behind systems consisting solely of rocky-icy bodies. In cases with r giant < r disk (upper- left 
region of Fig. 2) planets are born into a substantial gas disk, and such systems generally pro- 
duced a number of gas giants that migrated inward significantly. Planet-bearing systems can 
be further classified into "planet-dominated" ones, wherein the planets clear the disk from the 
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inside out, and "disk-dominated" cases, in which the disk clears the planets, their ultimate fate 
depending on what happens at the star-disk interface (not resolved in our simulations). The 
system depicted in Fig. 1 represents an intermediate case, with more examples shown in the 
Supporting Material. Whether planets accumulate at the original inner disk edge (|22ll23]|24l) 
or in an inner hole, migration tends to produce crowded systems, leading to the excitation of 
eccentricities via planet-planet scattering and resonances (I731 I21I31I7?!). We do not know where 
the true distribution of disk properties falls, but because 6-7% of Sun-like stars are observed to 
harbor giant planets (25), at least that fraction of disks must fall into the giant planet-bearing 
upper-left region of Fig. 2. Despite the speed of the hybrid code, it is as yet not computa- 
tionally feasible to perform simulations of a sufficient fidelity and number to allow a detailed 
comparison to the statistical properties of the discovered exoplanets. Instead, we performed a 
smaller set of lengthier simulations, focusing principally on a part of the planet-bearing region. 
These produced qualitative agreement with some of the key features of the exoplanet ensemble, 
namely the correlation of host-star metallicity with planet occurrence (|26]|27l), the distribution 
of planet masses, planet orbital periods, and the mass-eccentricity distribution (|28|) . as well as 
insights into how these features arose (see Supporting Material for details). 

Our results also suggest how the Solar System fits into the picture. In systems with r gian t ~ 
Tdi S k, gas giants do form, but undergo only modest migration and eccentricity growth; thus, it 
is here where we would most naturally expect to find a Solar System-like outcome. Figure 2 
shows that these cases occupy a relatively narrow region within the parameter space, roughly 
a diagonal line extending from (a = 1 x 10~ 3 , Mdi S k = 0.03Ms UI1 ) to (a = 10~ 2 , Mdi S k = 
0.08M Sun ). Thus, whatever the true distribution of disks within Fig. 2 — unless it just happens 
that disks with r g i an t ~ Tdi S k are somehow preferred (29) — it is likely that only a minority will 
lie within this region. Furthermore, even within this subset there are large stochastic variations, 
as evidenced by Fig. 2; in only one of the outcomes (a = 3 x 1CT 3 , M disk = 0.05M Sun ) do 
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the gas giants bear a reasonable resemblance to Jupiter and Saturn (for an animated example 
see Supporting Material). The other potential pathway to a Solar System analogue are cases in 
which all gas giants except two are engulfed by the star; however, since both survivors must 
somehow themselves undergo little migration, such outcomes also appear improbable (Fig. 2 
shows no candidates for this scenario). All of this leads us to predict that within the diverse 
ensemble of planetary systems, ones resembling our own are the exception rather than the rule. 
Observations may be hinting at this already (30), though the true planet distribution remains 
largely obscured by selection effects (|25| . On the other hand, scaled-down versions of the Solar 
System, in which a moderate amount of migration took place, are likely to be more common; 
indeed, such a system has recently been discovered via microlensing 071) . Finally, scenarios 
in which type II migration is reduced d2I]|221) would modify our prediction, permitting a more 
common occurrence of Solar System analogues. 

In all of our simulations, the formation of a gas giant brings with it violent scattering of 
neighboring smaller bodies, including other cores about to undergo runaway gas accretion 
themselves. Such scattering has been proposed as the origin of Uranus and Neptune (34}), with 
dynamical friction from the remnant outer planetesimal disk (not modelled here) serving to pre- 
vent their ejection and ultimately re-circularize their orbits. Thus, whether or not Jupiter and 
Saturn analogues are rare, it is likely that Uranus and Neptune analogues are quite common. 
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Figure 1 : Example of a planet formation simulation. A Solar-metallicity disk with initial mass 
Mdisk = 0.088Ms un and viscosity parameter a = 7.5 x 1CT 3 , containing protoplanets of initial 
mass 1CT 3 Me between 2 and 30 AU, evolves for 0.5 Gyr. Top Left: Planet semi-major axis 
over time, with the azimuthally-averaged gas disk surface density overlaid as a contour plot. 
Right, top four panels: Snapshots of the system at different intermediate times, showing planet 
eccentricity (right scale) and inclination (far-right scale, indicated by a horizontal tickmark 
connected to each planet by a vertical line) versus semi-major axis. Planet cores (black empty 
circles) together with whatever gas envelopes they have accreted (dark blue) are labelled with 
their mass in Me (black and dark blue, respectively). Planets crossing the inner simulation 
boundary at 0.1 AU are removed, and their final mass and orbital elements shown (empty red 
circles). The azimuthally-averaged disk surface density is also shown (light blue, left scale). 
One planet ends up at the inner disk edge (likely < 0.1 AU if due to the star's magneto spheric 
cavity (|22l) : however, we set it at 0.2 AU for computational reasons). Bottom: The state of the 
system at 0.5 Gyrs. Orbital elements are averages over the last 1 Myr. 
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Figure 2: Final outcomes of a set of 100 simulations, spanning 1CT 3 to 1CT 2 in viscosity pa- 
rameter a, and 0.01Ms un to 0.1Ms un in initial disk mass. Simulations are ended after 0.5 Gyrs 
have elapsed, or if they fail to produce any gas giants within the gas disk's lifetime. Planet 
semi-major axes, eccentricities and inclinations (denoted as in Fig. 1) are plotted, as well as the 
relative solids and gas content of each planet (solid core: black; gaseous envelope: dark blue; 
size oc mass 1 / 3 , see mass scale at top). To keep computational cost reasonable, simulations 
have an inner boundary at 0.25 AU, beyond the initial inner edge of the gas disk; any body that 
crosses the boundary is removed, and a red circle is plotted showing its orbital elements and 
mass at the time of removal. Toward high M disk and low a, planets form early and often during 
the gas disk's lifetime, most migrate extensively, and many acquire high eccentricities in the 
process (See Fig. 1 and Supporting Material). Toward low M disk and high a, planet formation 
is too slow to produce any gas giants during the disk lifetime. Between these two extremes is a 
relatively narrow boundary region (thick borders) in which gas giants migrate little and remain 
at low eccentricity, thus producing some outcomes more similar to the Solar System. No gas 
giants at all form in disks of 0.02Ms un or less; for comparison, this is the approximate lower 
limit on the Solar System's birth disk, called the "minimum mass Solar nebula" (TJTj) . 



11 




1 10 100 

r (AU) 



Figure 3: Approximate timing and location of gas giant formation in a protoplanetary disk. 
Bottom panel: The final or "isolation" mass of solid cores (black dots), with spacing between 
successive cores taken from planet formation simulations diffl) . Top: The time (thick solid 
curve) for a core (black dots; vertical dotted lines connect to corresponding core in bottom 
panel) to become a gas giant (horizontal dotted lines show times for individual protoplanets). 
We approximate this as the sum of the time for the core to reach its final mass, r corc (thin solid 
curve), and the time for the core to undergo runaway gas accretion, taken to be its Kelvin- 
Helmholtz time (T39|) . t K h (dashed curve). As in more detailed calculations we find that 
gas giant formation commences at one particular radius which for typical parameters lies in or 
near the Jupiter- Saturn region; in this case at 7 AU, and at time r gian t just under 2 Myrs). Giant 
formation begins in a burst, with several planets growing in rapid succession, then slows down 
as it spreads to larger and smaller radii. In practice, once an inner hole forms in the gas disk, 
formation is constrained to progress only outwards) 
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